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Abstract 

In  this  work,  an  isothermal,  steady- state,  three-dimensional  (3D)  multicomponent  transport  model  is  developed  for  proton  exchange  membrane 
(PEM)  fuel  cell  with  straight  gas  channels.  The  model  computational  domain,  includes  anode  flow  channel,  membrane  electrode  assembly  (MEA) 
and  cathode  flow  channel.  The  catalyst  layer  within  the  domain  has  physical  volume  without  simplification.  A  comprehensive  set  of  3D  continuity 
equation,  momentum  equations  and  species  conservation  equations  are  formulated  to  describe  the  flow  and  species  transport  of  the  gas  mixture  in 
the  coupled  gas  channels  and  the  electrodes.  The  electrochemical  reaction  rate  is  modified  by  an  agglomerate  model  to  account  for  the  effect  of 
diffusion  resistance  through  catalyst  particle.  The  activation  overpotential  is  predicted  locally  in  the  catalyst  layer  by  separately  solving  electric 
potential  equations  of  membrane  phase  and  solid  phase.  The  model  is  validated  by  comparison  of  the  model  prediction  with  experimental  data  of 
Ticianelli  et  al.  The  results  indicate  the  detailed  distribution  characteristics  of  oxygen  concentration,  local  current  density  and  cathode  activation 
overpotential  at  different  current  densities.  The  distribution  patterns  are  relatively  uniform  at  low  average  current  density  and  are  severely  non- 
uniform  at  higher  current  density  due  to  the  mass  transfer  limitation.  The  local  effectiveness  factor  in  the  catalyst  layer  can  be  obtained  with  this 
model,  so  the  mass  transport  limitation  is  displayed  from  another  point  of  view. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  (PEM)  fuel  cells  are  supposed 
to  be  the  most  promising  candidate  for  powering  of  electric 
vehicles  due  to  their  high  power  density,  short  response  time, 
low  operating  temperature  and  pollution  free.  Their  scalability 
and  relatively  flexibility  in  terms  of  the  fuel  makes  them  prime 
candidates  for  a  variety  of  stationary  applications,  including 
fixed  power  generations,  distributed  power  systems  and  portable 
electronic  appliance.  Modeling  and  simulation  are  being  used 
extensively  in  researches  and  industrial  applications  to  gain  bet¬ 
ter  understanding  of  the  fundamental  processes  and  to  optimize 
fuel  cell  designs  before  building  a  prototype  for  engineering 
application. 

Over  the  past  few  years,  several  models  of  PEM  fuel  cells 
were  proposed  in  the  literature.  Bemardi  and  Verbrugge  [1,2] 
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and  Springer  [3]  proposed  one-dimensional  models  that  pro¬ 
vided  good  preliminary  foundations  for  PEM  fuel  cell  modeling. 
However,  a  one-dimensional  model  cannot  simulate  the  decrease 
of  reactants  and  the  accumulation  of  products  in  the  flow  direc¬ 
tion.  The  two-dimensional  models  by  Fuller  and  Newman  [4] 
and  Nguyen  and  White  [5]  assumed  that  diffusion  was  the  only 
mechanism  for  oxygen  transport  and  did  not  consider  the  inter¬ 
action  between  the  flow  in  the  channel  and  gas  diffusion  layer 
(GDL),  which  is  seemingly  a  drawback  and  the  model  needs  to 
be  further  refined.  More  recently,  a  general  trend  can  be  observed 
to  apply  the  methods  of  computational  fluid  dynamics  to  fuel  cell 
modeling.  By  using  computational  fluid  dynamics  (CFD),  Gurau 
and  Liu  [6]  presented  the  first  unified  approach  by  coupling  the 
flow  and  transport  governing  equations  in  the  gas  channel  and 
GDL.  In  their  model,  it  was  assumed  that  the  catalyst  layer  was 
infinitesimally  thin  and  the  process  in  the  catalyst  was  totally 
neglected.  In  the  work  by  Um  et  al.  [7],  the  catalyst  layer  was 
solved  with  the  assumption  that  the  solid  (or  electronic)  phase 
potential  is  uniform  across  the  catalyst  layer.  Wang  et  al.  [8] 
classified  four  regimes  of  water  transport  and  distribution  in  the 
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Nomenclature 

As  specific  reaction  area  of  the  catalyst  layer  (m-1) 
Ac h  gas  channel  cross-section  area  (m2) 

Am  geometrical  area  of  the  membrane  (m2) 
c  molar  concentration  (mol  m-3) 

D  mass  diffusivity  (m2  s-1) 

F  Faraday  constant  (96,487  C  mol- 1 ) 

H  height  (m),  Henry  constant  (dimensionless) 
i  volumetric  current  density  (A  m-3) 

I  average  current  density  (A  m-2) 

y'o  exchange  current  density  (A  m-2) 

k  reaction  rate  constant  (ms-1) 

K  permeability  of  electrode  (m2) 

L  length  (m) 

M  molecular  weight  (dimensionless) 

Mj  Thiele  modulus  (dimensionless) 

p  pressure  (Pa) 

R  universal  gas  constant  (8.314  5  J  mol-1  K-1) 

S  source  term  of  equations 

T  temperature  ( K) 

v  velocity  vector  (ms-1) 

V  electrical  potential  (V) 

v  coordinate  (m) 

y  coordinate  (m) 

z  coordinate  (m) 


Greek  letters 

a  transfer  coefficient  of  electrochemical  reaction 

(dimensionless) 

/3  net  water  transport  coefficient  per  proton 

(dimensionless) 

s  porosity  of  electrode  (dimensionless) 

f  stoichiometric  flow  ratio  (dimensionless) 

0  viscosity  (kg  m- 1  s  - 1 ) ;  overpotential  (V) 

p  density  (kg  m- 3 ) 

a  electrical  conductivity  (S  m- 1 ) 

0  electrical  potential  (V) 

co  species  mass  fraction  (dimensionless) 

Subscripts 
av  average 

a  anode 

c  cathode 

ct  catalyst 

ch  channel 

eff  effective 

g  gas 

i  species 

in  inlet 

k  anode  or  cathode 

m  membrane 

max  maximum  momentum 

o  oxygen 

oc  open  circuit 

ref  reference  values 


s  solid;  speific 

tot  total 

w  water 

Superscripts 
a  anode 

c  cathode 

m  membrane 


PEMFC  air  cathode  and  presented  some  interesting  two-phase 
flow  and  transport  results.  Recently  a  comprehensive  3D  model 
of  PEM  fuel  cell  was  presented  by  Berning  et  al.  [9].  Again,  in 
their  model  the  catalyst  layer  was  treated  as  a  surface  without 
thickness.  The  same  assumption  had  been  made  in  the  modeling 
approach  presented  by  Dutta  et  al.  [10],  who  obtained  results  by 
adding  the  necessary  sources  and  correction  terms  to  the  govern¬ 
ing  equations  within  the  framework  of  commercial  CFD  codes. 
Other  notable  work  in  this  area  includes  models  developed  by 
N guy en  and  co-workers  [  1 1 , 1 2] .  In  their  models ,  the  catalyst  lay¬ 
ers  were  also  assumed  to  be  infinitesimally  thin.  Thus  it  can  be 
observed  that  because  of  the  very  complicated  physical  phenom¬ 
ena  occurring  in  the  PEM  fuel  cell  process,  in  the  existing  models 
proposed  so  far  different  assumptions  are  often  to  be  made  in 
order  to  make  the  problem  solvable.  The  three  major  assump¬ 
tions  are  related  to  the  following  aspects:  (1)  whether  the  anode 
or  cathode  are  both  included;  (2)  whether  the  activation  overpo¬ 
tential  is  assumed  to  be  a  constant,  and  (3)  whether  the  diffusion 
resistance  through  the  catalyst  particle  is  taken  into  account. 

In  this  paper,  the  assumptions  mentioned  above  are  discarded, 
and  a  more  comprehensive  model  is  proposed  and  numerically 
solved  with  a  code  developed  by  the  authors.  In  the  follow¬ 
ing  presentation,  the  model  description  will  first  be  stated  in 
detail,  including  the  governing  equations  and  the  related  bound¬ 
ary  conditions,  followed  by  a  brief  presentation  of  the  numerical 
algorithm  adopted,  then  detailed  discussion  will  be  made  on  the 
numerical  results.  Finally,  some  conclusions  will  be  drawn. 

2.  Model  description 

Numerical  simulation  is  made  for  the  PEM  fuel  cell  with 
straight  and  parallel  channels  of  the  polar  plate,  shown  in  Fig.  1 . 
The  multi-channel  structure  and  the  periodicity  of  the  polar  plate 
makes  it  possible  to  select  one  unit  (shown  by  the  dashed  lines  in 
Fig.  1)  as  the  representative  one  of  the  whole  fuel  cell.  It  includes 
anode  gas  channel,  anode  GDL,  anode  catalyst  layer,  membrane, 
cathode  catalyst  layer,  cathode  GDL  and  cathode  gas  channel. 
Again  because  of  the  symmetry  of  the  left  and  right  half  unit, 
we  take  the  right  half  unit  as  our  computational  domain.  For 
the  clarity  of  presentation,  the  two-dimensional  cross-sections 
of  y-z  plane  and  x-z  plane  are  shown  in  Fig.  2.  Both  the  left  and 
right  boundaries  (in  y-direction)  of  Fig.  2a  are  the  symmetric 
surface. 

Reactants  move  from  the  gas  channel  into  the  GDL,  which 
serves  to  make  more  uniform  distribution  of  the  reactants  over 
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Fig.  1.  Straight  and  parallel  flow  field  and  flow  channel  (the  shaded  area  is  the 
polar  plate). 


the  catalyst  layer.  In  the  catalyst  layer,  the  reactants  are  trans¬ 
ported  by  diffusion  and  advection  to  participate  in  the  electro¬ 
chemical  reaction.  The  membrane  can  transport  the  proton  and 
dissolved  water  but  it  is  assumed  to  be  impermeable  for  gas. 

2.7.  Model  assumptions 

In  order  to  make  the  numerical  simulation  manageable,  some 
assumptions  are  made  as  follows: 

(1)  The  cell  operates  under  steady-state  condition. 

(2)  The  cell  temperature  is  uniform  and  fixed. 

(3)  The  reactant  and  production  are  assumed  to  be  ideal  gas 
mixtures.  And  the  produced  water  is  treated  as  vapor. 

(4)  The  electrode  is  treated  as  an  isotropic  and  homogenous 
porous  medium  and  the  properties,  such  as  porosity  and 
permeability  are  constants. 

(5)  The  membrane  is  impermeable  for  gas  phase. 

(6)  Ohmic  losses  in  the  GDL  and  current  collector  (or  bipolar 
plate)  are  neglected. 

(7)  The  flow  in  the  channels  is  considered  laminar. 

2.2.  Governing  equations 


(a)  y-z  cross  section 


O2,  N2 


HoO 


O2  N2 


(b) 


x-z  cross  section 


ssr  a 

At 


Fig.  2.  The  two-dimensional  cross-sections  of  the  computational  domain:  (a) 
y-z  cross-section  and  (b)  x-z  cross-section  ((1)  anode  gas  channel,  (2)  anode 
diffusion  layer,  (3)  catalyst  layer,  (4)  membrane,  (5)  cathode  diffusion  layer,  (6) 
cathode  gas  channel,  (7)  polar  plate). 


The  transport  of  gas  mixtures  in  the  gas  channels  and  in  the 
electrodes  conforms  to  the  mass,  momentum,  and  species  con¬ 
servation  principles.  The  corresponding  governing  equations  are 


written  as  follows: 

•  Mass  conservation  equation 

V  •  (pu)  =  (1) 

•  Momentum  conservation  equation 

1/eV  •  (pirn)  =  — Vp  +  1/eV  •  (77V u)  —  (rj/k) u  (2) 

•  Species  mass  fraction  conservation  equations 

V  •  (puu%)  =  V  •  (pAi,effVu%)  +  Sh  (3) 

V  •  (pu<u0)  =  V  •  (p7)o  effV(Uo)  +  S0  (4) 

V  •  (puo%)  =  V  •  (p7)w  effVcuw)  +  Sw  (5) 


where  u  is  the  superficial  velocity  vector,  which  is  propor¬ 
tional  to  the  fluid  velocity  vector  by  a  coefficient  e,  i.e.  poros¬ 
ity  of  the  porous  electrodes,  coi  the  species  mass  fraction  of 
the  gas  mixture,  A?eff  the  effective  species  diffusivity,  and  Si 
is  the  source  term  of  equations.  In  a  pure  fluid  region,  e  is 
unity  and  then  the  superficial  velocity  vector  is  reduced  to  the 
real  fluid  velocity  vector. 

The  above  governing  equations  are  assumed  to  be  applica¬ 
ble  for  both  gas  channels  and  electrodes.  Thus,  the  interfacial 
conditions  at  the  interfaces  of  channel-gas  diffusion  layer,  gas 
diffusion-catalyst  layer  and  catalyst-membrane  layer,  are  not 
needed.  Some  explanations  are  given  below  for  the  implemen¬ 
tation  of  above  governing  equations  in  individual  part  of  the 
entire  computational  domain.  In  the  general  form  of  momentum 
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conservation  equation,  Eq.  (2)  the  last  term  on  the  right  side  rep¬ 
resents  Darcy’s  drag  force  imposed  by  the  pore  walls  on  the  fluid 
within  the  pores,  which  usually  results  in  a  significant  pressure 
drop  across  the  porous  medium.  It  is  often  called  as  the  micro 
scale  viscous  term  or  Darcy’s  viscous  term.  Thus  in  the  porous 
medium,  the  general  momentum  conservation  equation  reduces 
to  the  extended  Darcy’s  law  for  the  flow  in  porous  media  with 
a  small  permeability.  While  inside  a  pure  fluid  region,  i.e.  gas 
channel,  it  recovers  the  standard  Navier-Stokes  equation  with 
the  porosity  being  unity  and  the  permeability  being  infiniteness. 
The  last  terms  in  Eq.  (1)  and  (3)— (5)  are  the  volumetric  sink  or 
source  terms  due  to  the  electrochemical  reactions  in  the  cata¬ 
lyst  layer,  and  they  are  zero  in  other  parts  of  the  computational 
domain. 

The  thermophysical  properties  of  the  mixtures  in  the  above 
governing  equations  are  determined  as  follows: 

•  The  density  of  gas  mixture  is  estimated  by  the  ideal  gas  law 

p  =  Mp/RT  (6) 

M=  (7) 

i 

where  M  and  M;  are  the  molecular  weights  of  gas  mixture  and 
species,  respectively,  and  R  is  the  universal  gas  constant. 

•  The  viscosity  of  the  mixture  is  specified  by 

ii  =  y^m  (8) 

i 

where  rj  and  rji  are  the  viscosities  of  the  gas  mixture  and 
species,  respectively. 

•  The  diffusivity  of  the  species  is  determined  for  flow  in  the 
porous  media  by  the  so-called  Bruggeman  model  [9] 

A,e  ff  =  Die1-5  (9) 

where  D*  is  the  diffusivity  of  gas  species  in  a  nonporous  sys¬ 
tem.  According  to  [9],  it  is  related  with  the  diffusivity  under 
the  reference  condition  expressed  by 

Di  =  Dljef(T/TKf)L5(p/prW~'  (10) 

The  non-zero  source  terms  existing  in  Eqs.  (1)  and  (3)— (5) 
for  the  two  catalyst  layers  are  given  by: 


•  Cathode  catalyst  layer 


So  =  -0c/4F)Mo 

(11) 

Sw  =  [(l  +2p)i0/2F]Mw 

(12) 

Sm  =  S0  +  Sw 

(13) 

Anode  catalyst  layer 

Sh  =  -(/a/2  F)Mh 

(14) 

Sw  =  -(fiU/F)Mv 

(15) 

Sm  —  Sh  +  sw 

(16) 

where  i  is  the  local  current  density  in  the  catalyst  layer,  F  the 
Faraday  constant,  and  /3  is  the  net  water  transport  rate  through 
the  membrane  per  proton. 

The  local  current  densities  in  anode  and  cathode  can  be 
obtained  by  the  Bulter-Volmer  equation.  Due  to  the  different 
characteristics  of  anodic  and  cathodic  electrochemical  reaction, 
the  polarization  potential  loss  is  minor  in  anode  and  relatively 
great  in  cathode.  Also  the  effect  of  the  reactant  concentration 
on  the  reactive  rate,  i.e.  current  density,  should  be  taken  into 
account.  So  the  original  Bulter-Volmer  equation  is  modified 
for  calculating  the  anodic  and  cathodic  local  current  density, 
expressed  as  follows: 

ia  =  As7o«/Cef)1/2  '  C n^FriJRT )  (17) 

k  =  Asj%(c™  /c^ref)  exp (ancFric/RT)  (18) 

where  jQ  is  the  exchange  current  density,  r)a  and  0c  are  activation 
overpotentials  of  the  anode  and  cathode,  respectively,  cJJ1  and  c™ 
are  the  molar  concentrations  of  hydrogen  and  oxygen  dissolved 
in  the  membrane  (indicated  by  the  superscript  m)  phase,  respec¬ 
tively,  a  the  cathode  transfer  coefficient,  na  the  electron  number 
of  anode  reaction  and  nc  is  that  of  cathode  reaction. 

The  dissolved  molar  concentration  of  species  in  the  polymer 
phase  is  given  by  Henry’s  law  [14]: 

cf1  =  Hc{  (19) 

where  H  is  the  Henry  constant,  q  and  c f1  are  the  molar  concen¬ 
trations  of  species  existing  in  gas  phase  and  membrane  phase  (or 
Nation,  polymer  phase),  respectively,  and  the  subscript  i  denotes 
species,  such  as  hydrogen  or  oxygen. 

Attention  is  now  turned  to  the  diffusion  process  in  the  cat¬ 
alyst.  Research  has  shown  that  the  catalyst  layer  is  porous  and 
made  up  of  clumps  of  carbon- supported  Pt  catalyst,  surrounded 
by  a  thin  layer  of  Nation  [14-16].  These  clumps  are  referred  to 
as  agglomerates.  The  gaseous  reactants  must  dissolve  into  the 
polymer  phase  and  diffuse  through  the  polymer  film  to  reach 
the  reaction  sites.  In  order  to  account  for  the  effect  of  diffu¬ 
sion  resistance  through  the  catalyst  with  porous  and  agglomerate 
structure,  the  local  current  density  is  modified  by  an  effective¬ 
ness  factor,  0 ,  which  is  a  measure  of  how  readily  reactants  diffuse 
through  the  catalyst  particle: 

i  =  0  •  i  (20a) 

According  to  [17],  0  can  be  determined  by  following  equation: 

6  =  tanh  Mj/Mj,  MT  =  Lcl  Jk/ Of  (20b) 

where  Mj  is  called  Thiele  modulus,  Lct  the  characteristic  length 
of  catalyst  particle,  k  the  reaction  rate  constant  and  Df1  is  species 
diffusivity  of  reactant  in  the  polymer  phase. 

The  reaction  rate  constant  can  be  expressed  as 

k  •  (cff  =  i/nkF  (21) 

where  n  is  the  order  of  reaction,  which  is  unity  for  cathode  reac¬ 
tion  and  1/2  for  anode  reaction,  and  nk  is  the  number  of  electron 
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transferred  in  the  anodic/cathodic  electrochemical  reaction  (sub¬ 
script  k  denotes  anode  or  cathode). 

An  effectiveness  factor  of  1.0  indicates  that  reactants  diffuse 
through  the  agglomerate  catalyst  particle  without  resistance.  The 
factor  less  than  1.0  represents  that  the  agglomerate  offers  some 
resistance  to  reactant  diffusion,  thereby,  limiting  the  reaction 
rate.  This  corrected  method  is  so-called  the  agglomerate  model 
in  literatures  [14-16].  During  the  computation,  the  reaction  rate 
constant  k  is  calculated  from  the  local  current  density  and  molar 
concentration  according  to  Eq.  (21). 

The  above  descriptions  are  made  for  the  formulation  of  the 
gas  mixture  transport  in  gas  channels  and  porous  electrodes. 
Attention  is  now  turned  to  the  electrical  potential  prediction. 
For  the  MEA,  two  potential  equations  are  solved  in  the  present 
model.  The  solid  phase  potential  equation  represents  transport 
of  electrons  in  the  solid  conductive  regions  (GDL  and  catalyst 
layer),  which  reads: 


V  •  (<rsV0s)  +  S<p,s  =  0  (22) 

The  membrane  phase  potential  equation  depicts  transport  of  pro¬ 
tons  in  the  MEA  that  consists  of  both  catalyst  layers  and  the 
membrane  itself,  expressed  by: 

V  ’  (crmV0m)  T  S(p,m  —  0  (23) 

In  the  above  equations  0S  and  0m  denote  electrical  potential 
of  solid  phase  and  membrane  phase  respectively,  as  the  electrical 
conductivity  of  the  solid  phase  and  crm  is  the  protonic  conduc¬ 
tivity  of  the  membrane  phase.  The  terms  S^s  and  S^m  are  the 
volumetric  source  terms,  which  exist  only  in  the  catalyst  layer 
and  are  determined  based  upon  the  transfer  current  densities  as 
follows: 


•  For  cathode  catalyst  layer 


—  A 

(24) 

Sf,  s  —  ic 

(25) 

For  anode  catalyst  layer 

absolute  value  of  the  difference  between  the  reference  potential 
and  the  predicted  potential  at  that  location  [18]. 


Oele  — 


Opro  — 


0  s  0s, ref 

0m  0m, ref 


(28) 

(29) 


where  rje le  is  electric  Ohmic  overpotential  and  f/pro  is  protonic 
Ohmic  overpotential. 

The  activation  overpotential  is  then  obtained  by: 


Osl  —  7?a,tot  ?7a,ele  ^a,pro  (30) 

0c  =  0c,tot  —  0c,clc  —  7?c,pro  (31) 

where  rj &,tot  and  rjc, tot  are  the  total  overpotential  including  all 
potential  losses  of  anode  and  cathode,  respectively.  In  this  model, 
rj Cjtot  is  assumed  to  be  a  known  quantity  and  rj &,tot  is  obtained 
in  the  calculation  to  ensure  the  anode  average  current  density 
equal  to  that  of  the  cathode. 

The  average  current  density  is  computed  by  integrating  the 
local  current  densities  over  all  control  volumes  in  the  catalyst 
layer.  It  is  necessary  to  make  the  anode  average  current  density 
equal  to  that  of  cathode  because  of  conservation  of  the  electric 
current,  i.e.: 

I  =  (1  /Am)  J2(k  ■  VCv)  =  (1  /Am)  5>  •  Vcv)  (32) 

where  Am  is  the  geometrical  area  of  the  membrane  and  Vcv  is 
the  volume  of  unit  control  volume  grid. 

The  operating  potential  of  the  cell  is  then  calculated  by: 


Ecell  —  Vqc  ?7a,tot  7?c,tot  7?m,pro  (33) 

where  Voc  is  the  open  circuit  potential  of  the  cell  and  f/m,pro  is 
the  Ohmic  overpotential  in  the  membrane. 

The  open  circuit  potential  is  usually  obtained  by  Nernst  equa¬ 
tion.  Examination  of  Nernst  equation  shows  a  decrease  of  the 
open  circuit  potential  with  temperature.  But  experimental  results 
of  Parthasarathy  et  al.  [19]  show  an  opposite  effect.  Thus  we  use 
the  empirical  results  of  Parthasarathy  et  al.,  which  were  fitted 
by  a  linear  function  of  temperature  in  [6]. 


(26)  Voc  =  0.025  T  +  0.2329 


(34) 


(27) 


2.3.  Boundary  conditions 


In  order  to  be  able  to  calculate  the  overpotential  in  the  dif¬ 
ferent  regions,  it  is  important  to  select  a  reference  point  where 
the  electrical  potential  is  zero.  As  usual  we  define  the  electri¬ 
cal  potential  at  the  interface  between  the  bipolar  plate  and  the 
anode  GDL  to  be  zero,  which  is  called  the  reference  potential 
hereafter.  The  fuel  cell  operating  potential  is  then  the  potential  at 
the  interface  between  the  bipolar  plate  and  the  cathode  GDL.  The 
geometrical  interface  between  the  anode  catalyst  layer  and  the 
membrane  is  selected  as  the  reference  plane  for  the  membrane 
phase  potential. 

We  then  give  the  definition  of  the  Ohmic  type  over  potential 
(or  Ohmic  loss),  including  electric  and  protonic  Ohmic  over¬ 
potential.  The  over  potential  in  a  local  place  is  defined  as  the 


At  the  y-z  plane  boundaries,  symmetric  is  assumed,  i.e.  all 
gradients  in  the  y-direction  are  set  to  zero  at  these  two  boundary 
planes  of  the  domain  (Fig.  2a). 

At  the  x-y  planes  boundaries,  i.e.  the  interfaces  between  the 
gas  channel  and  current  collector  plates  (Fig.  2b),  the  no-slip  and 
impermeability  conditions  are  implemented  in  the  z-direction. 

The  boundary  conditions  at  the  two  y-z  planes  are  now  dis¬ 
cussed. 

At  the  gas  channel  inlet,  the  gas  mixture  average  velocity 
and  the  species  concentration  are  prescribed.  The  inlet  velocity 
is  specified  by: 

^k,in  —  £k(7max/  ^kE)(7?7in/ /?k,in)(l/  ^k,in)(Am/ Ach) 


(35) 
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where  is  the  reactant  stoichiometric  flow  ratio  of  anode  or 
cathode,  subscript  k  denotes  the  electrode,  7max  the  maximum 
average  current  density,  Am  the  geometrical  area  in  x-y  plane  of 
the  membrane,  and  Ac h  is  the  cross-section  area  of  gas  channel 
in  y-z  plane. 

At  the  gas  channel  outlet,  the  pressure  is  prescribed  as 
the  operating  pressure.  For  all  other  dependent  variables,  their 
change  rates  are  assumed  infinitesimal,  i.e.  the  gradients  in  the 
v-direction  are  set  to  zero. 

The  gas  diffusion  layer  and  catalyst  layer  are  surrounded 
by  the  sealed  plates  at  inlet  and  outlet  plane.  So  the  bound¬ 
ary  conditions  of  gas  diffusion  layer  and  catalyst  layer  at  inlet 
and  outlet  planes  take  non- slip  condition  for  the  velocity  and 
non-permeable  condition  for  the  species  mass  fraction. 

There  is  no  need  to  give  the  boundary  conditions  at  the  inter¬ 
face  between  different  parts  of  cell,  such  as  the  interface  between 
channel  and  GDL,  and  the  interface  between  the  GDL  and  the 
catalyst  layer,  since  they  are  within  the  computational  domain. 
And  it  is  in  this  sense  that  makes  our  computation  conjugated. 
The  multi  domain  technique,  which  includes  two  gas  channels, 
two  catalyst  layers,  two  gas  diffusion  layers  and  membrane,  is 
used  in  the  numerical  simulation.  So  the  membrane  is  treated  as 
the  fluid  with  infinite  viscosity,  which  can  ensure  that  the  velocity 
in  the  membrane  equals  zero  [20] .  For  the  species  concentration 
equation,  the  species  mass  diffusivity  in  the  membrane  is  set  to 
zero.  Thus  non-slip  condition  for  the  velocity  and  non-permeable 
condition  for  the  species  concentration  are  implemented  in  the 
numerical  simulation  for  the  interface  between  membrane  and 
catalyst  layer. 

3.  Numerical  algorithm  and  methods 

The  governing  equations  (Eqs.  (1)— (5),  (22),  (23)),  together 
with  the  boundary  conditions  are  discretized  by  the  finite- 
volume  method.  The  SIMPLE  algorithm  of  Patankar  and  Spald¬ 
ing  [20,21]  is  utilized  to  deal  with  the  coupling  of  the  velocity 
and  field.  Since  all  governing  equations  are  coupled  with  each 
other  through  the  source  terms,  they  ought  to  be  solved  simulta¬ 
neously  with  the  iterative  method.  The  solution  procedure  is  as 
follows.  At  first,  the  initialization  of  the  dependent  variables  is 
assumed.  The  initial  velocity  and  concentration  (or  mass  frac¬ 
tion)  take  the  inlet  values  for  the  gas  channel  and  a  small  value 
in  the  electrode,  such  as  10-3.  The  potential  is  initialized  as 
zero  or  a  small  value  for  whole  computational  domain  (includ¬ 
ing  the  membrane  and  two  catalyst  layers).  Then  the  velocity 
and  pressure  fields  for  the  gas  mixture  are  solved  in  the  cou¬ 
pled  gas  channel  and  porous  electrode  domains  with  previous 
(or  assumed)  current  densities,  followed  by  the  solution  of  the 
gas  species  mass  fraction  equations,  which  are  also  dependent 
on  the  local  current  densities.  Then  the  potential  governing  equa¬ 
tions,  Eqs.  (22)  and  (23),  are  solved  for  calculating  the  Ohmic 
overpotential,  and  the  activation  overpotential  can  be  obtained 
by  Eqs.  (30)  and  (31).  After  this,  the  local  current  density  is 
computed  from  the  solved  values  of  reactant  species  concentra¬ 
tion  and  activation  overpotential  according  to  the  Bulter-Volmer 
equation.  The  newly  solved  current  densities  and  other  variables 
are  compared  with  the  previous  ones.  If  the  convergence  con- 


Table  1 


The  design  and  operating  parameters  of  the  cell 


Parameters 

Value 

Gas  channel  length  (L) 

0.07112m 

Gas  channel  width  (W) 

7.62  x  10-4  m 

Gas  channel  height  (Hc h) 

7.62  x  10-4  m 

Diffuser  layer  height  (Hf) 

2.54  x  10“4  m 

Catalyst  layer  height  ( Hct ) 

2.87  x  10~5  m 

Membrane  height  (Hm) 

2.3  x  10~4  m 

Temperature  ( T) 

353  K  [22] 

Anode/cathode  pressure  ( pjpc ) 

3/5  atm.  [22] 

Fuel  stoichiometric  flow  ratio  (£a) 

3 

Air  stoichiometric  flow  ratio  (£c) 

3 

Relative  humidity  of  inlet  fuel  (RHa) 

100% 

Relative  humidity  of  inlet  air  (RHC) 

0 

Oxygen  mass  fraction  of  inlet  air  (co0) 

0.232 

Gas  diffusion  layer  porosity  (Ai) 

0.4  [2] 

Catalyst  layer  porosity  (sc) 

0.28 

Characteristic  length  of  catalyst  (Lct) 

1.0  x  10"6m  [14] 

dition  is  not  satisfied,  the  above  solution  procedure  is  repeated 
(next  level  of  iterations). 

The  solution  is  considered  to  be  convergent  when  the  rela¬ 
tive  error  of  each  dependent  variable  between  two  consecutive 
iterations  is  less  than  10~5. 

In  conducting  a  numerical  simulation,  a  great  number  of 
parameters  and  physical  properties  are  required.  They  are  all 
listed  in  Tables  1  and  2.  These  values  are  basically  adopted 
from  [22]  whose  experimental  data  are  used  to  compare  with 
our  numerical  prediction.  Part  of  the  data  which  were  not  sup¬ 
plied  in  [22]  are  adopted  from  some  similar  Refs.  [6-10].  In 

Table  2 

The  physical  properties  used  in  the  model 

Properties 

Value 

Permeability  of  diffuser  ( K d) 

1.76  x  10“ 11  m2  [6] 

Permeability  of  catalyst  layer  ( Kc ) 

1  x  10“14m2 

Reference  diffusivity  of  H2  in  gas 

0.915  cm2  s_1  (1  atm,  307  K)  [9] 

(Dh,ref) 

Reference  diffusivity  of  O2  in  gas 

0.220  cm2  s_1  (1  atm,  293  K)  [9] 

(Do, ref) 

Reference  diffusivity  of  H2O  in  gas 

0.256  cm2  s_1  (1  atm,  307  K)  [9] 

(Dw,ref) 

Henry  constant  of  H2  in  the  Nation 

0.19  [14] 

(H) 

Henry  constant  of  O2  in  the  Nation 

0.64  [14] 

(H) 

H2  reference  concentration  in  the 

56.4  mol  m-3  [2] 

Nafion  (c”ref) 

O2  reference  concentration  in  the 

3.39  mol m-3  [2] 

Nafion  (c“ef) 

Diffusivity  of  O2  in  the  Nafion  (D™ ) 

1.22  x  10“ 10  m2  s-1  [2] 

Diffusivity  of  H2  in  the  Nafion  (DjJ1) 

2.59  x  10“ 10  m2  s-1  [2] 

Cathode  transfer  coefficient  (a) 

0.5 

Exchange  current  density  multiply 

2.0  x  108  Am-3 

specific  area  for  anode  (As  j'q) 

Exchange  current  density  multiply 

1.6  x  102  Am-3 

specific  area  for  cathode  (As  jcf) 

Solid  phase  conductivity  (crs) 

53  Sm-1  [2] 

Membrane  phase  conductivity  (crm) 

17  S  m-1  [2] 

Net  water  transfer  rate  (/ 3 ) 

0.2  [5] 
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Table  3 


Grid  independence  test 


Stage 

Grid  size 

4v  (A  cm  2) 

Percentage  change  (abs.) 

1 

30  x  20  x  36 

0.39641 

— 

40  x  20  x  36 

0.39642 

0.0025 

50  x  20  x  36 

0.39643 

0.0025 

2 

40  x  10  x  36 

0.39621 

— 

40  x  20  x  36 

0.39642 

0.053 

40  x  30  x  36 

0.39654 

0.03 

3 

40  x  20  x  28 

0.39637 

— 

40  x  20  x  36 

0.39642 

0.013 

40  x  20  x  46 

0.39646 

0.01 

addition,  the  total  cathodic  overpotential  is  assumed  as  a  known 
quantity  and  then  the  average  current  density  and  cell  operat¬ 
ing  potential  can  be  calculated  when  the  solution  of  governing 
equations  is  convergent.  Different  values  of  the  total  cathodic 
overpotential  are  applied,  so  the  polarization  curve  of  the  cell 
can  be  obtained. 

The  grid  system  used  is  40  (length)  x  20  (width)  x  36 
(height).  The  gas  channel  is  divided  into  eight  control  volumes 
in  the  z-direction,  the  GDL  5,  catalyst  layer  3  and  membrane 
4.  The  grid  independence  test  is  performed  by  increasing  and 
decreasing  the  number  of  the  grid  cells.  The  test  is  conducted 
in  three  stages:  (1)  grids  of  y-  and  z-directions  are  fixed  while 
grid  of  v-direction  is  varied,  (2)  grids  of  x-  and  z-directions  are 
fixed  while  grid  of  y-direction  is  varied  and  (3)  grids  of  x-  and 
y-directions  are  fixed  while  grid  of  z-direction  is  varied.  The 
results  are  summarized  in  Table  3.  Considering  both  accuracy 
and  economics,  it  can  be  found  from  Table  3  that  the  grid  system 
of  40  x  20  x  36  is  appropriate  for  the  present  study.  Hence,  all 
the  subsequent  calculations  in  the  present  study  are  performed 
using  the  grid  system. 

4.  Results  and  discussion 

In  this  section,  the  predicted  V-I  curve  will  first  be  compared 
with  test  data  available  to  us,  and  the  possible  reasons  which  may 
account  for  some  discrepancy  be  discussed.  Then  the  predicted 
results  of  the  pressure  field,  the  gas  mixture  velocity  fields,  and 
the  distribution  of  oxygen  mass  fraction  in  the  cathode,  local 
current  density,  local  activation  overpotential,  and  local  effec¬ 
tiveness  factor  in  the  catalyst  layer  will  be  presented  in  order. 
Discussion  of  the  major  features  of  the  present  study  will  also 
be  provided. 

The  predicted  fuel  cell  polarization  curve  of  the  fuel  cell  stud¬ 
ied  is  shown  in  Fig.  3.  Provided  there  are  also  the  measurement 
results  of  Ticianelli  et  al.  [22],  from  which  the  major  operational 
parameters  are  adopted.  It  is  a  common  practice  in  literatures  that 
the  test  result  of  [22]  is  taken  as  a  kind  of  the  benchmark  data  for 
the  validation  of  a  numerical  model  [6-10].  The  predicted  curve 
agrees  with  the  measured  one  very  well  at  low  current  densities, 
when  the  performance  of  the  cell  is  governed  by  the  electrode 
kinetics,  and  most  of  potential  losses  are  due  to  electrode  acti¬ 
vation.  While  at  the  high  current  density  (>0.8  A  cm-2),  there 
is  some  discrepancy  between  the  predicted  results  and  experi- 


Fig.  3.  Model  prediction  of  the  polarization  cure  compared  with  experimental 
data. 

mental  data,  with  the  predicted  value  being  higher  than  that  of 
measured.  It  is  probably  due  to  the  model  assumption  that  there 
is  no  liquid  formed  in  the  electrode.  At  high  current  densities, 
a  significant  amount  of  water  is  produced  at  the  cathode.  Also 
more  water  is  transported  to  the  cathode  by  the  electro-osmotic 
force  from  the  anode.  The  excess  water  floods  the  cathode,  clogs 
the  pores  and  prevents  oxygen  transport  to  the  catalyst  layer, 
leading  to  the  deterioration  of  the  cell  performance.  Thus  in  the 
high  current  density  region,  the  present  model  overestimates  to 
some  extent  the  cell  performance,  oxygen  concentration  and  the 
local  current  density,  with  the  neglect  of  the  liquid  water  exist¬ 
ing  in  the  cathode.  However,  the  predicted  distribution  patterns 
of  the  dependent  variables  can  still  provide  useful  information 
for  further  understanding  the  complicated  process.  The  devel¬ 
opment  of  a  two-phase  water  transport  model  is  now  underway 
and  the  analysis  of  influence  of  liquid  water  flooding  on  cell 
performance  will  be  analyzed  in  our  future  work. 

Another  possible  reason  which  may  account  for  some  dis¬ 
crepancy  between  predicted  and  test  results  is  the  periodicity 
assumption.  The  test  data  was  obtained  for  an  entire  fuel  cell. 
The  real  fuel  cell  current  collecting  plate  is  machined  with  a 
finite  number  of  parallel  channels.  In  this  study,  one  half  typical 
unit  of  the  cell  is  selected  as  the  computational  domain  to  reduce 
the  computational  cost,  and  the  symmetric  boundary  condition 
is  adopted  in  the  y-direction.  Although  this  is  a  common  prac¬ 
tice  in  fuel  cell  simulation,  this  treatment  implies  that  the  fuel 
cell  is  infinite  in  y-direction.  Thus,  it  may  lead  to  some  error  in 
predicting  the  cell  performance. 

A  question  may  arise  as  whether  the  isothermal  model 
adopted  in  this  study  may  lead  to  some  discrepancy.  According 
the  simulation  results  in  [9,10],  the  temperature  rise  predicted 
in  the  cell  was  only  in  the  order  of  several  degrees.  Thus,  it  is 
expected  that  the  isothermal  model  is  appropriate  for  prediction 
of  the  cell  performance. 

Fig.  4  shows  the  relative  pressure  field  compared  to  the  outlet 
pressure  in  both  electrodes.  The  pressure  gradient  is  small  along 
the  channel  and  the  pressure  difference  between  the  inlet  and 
outlet  is  less  than  20  Pa.  There  is  a  small  pressure  drop  in  the 
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Fig.  4.  Relative  pressure  (p ,  Pa)  fields  in  the  electrodes:  (a)  anode  and  (b)  cath¬ 
ode. 

anode  while  a  small  raise  in  the  cathode.  It  is  due  to  the  mass 
sink  or  source  resulting  from  the  electrochemical  reaction  in  the 
catalyst  layer. 

The  gas  mixture  velocity  fields  in  the  anode  and  cathode  are 
shown  in  Figs.  5  and  6.  In  the  figure,  the  velocity  vectors  indicate 
both  the  direction  of  the  flow  and  the  velocity  magnitude  with  a 
scaling  velocity  adhering  in  each  plot.  Figs.  5a  and  6a  represent 
the  flow  in  the  y-z  cross-section  adjacent  to  the  middle  of  the 
channel  and  Figs.  5b  and  6b  the  flow  in  the  x-y  cross-section 
in  the  middle  of  the  diffusion  layer.  From  the  figures,  following 
features  may  be  noted.  Firstly,  it  can  be  clearly  observed  that 
the  gas  mixtures  move  much  more  slowly  in  the  two  porous 
electrodes  than  that  in  the  gas  channels.  Secondly,  the  scaling 
velocities  are  quite  different  for  the  cathode  and  anode  diffusion 
layers.  It  can  be  found  that  the  magnitude  of  velocity  vectors 
in  the  anode  diffusion  layer  is  one  order  larger  than  that  of  the 
cathode  one. 

The  Peclet  number  is  calculated  for  the  gas  flow  towards 
catalyst,  being  0.0058  for  cathode  side  and  0.0009  for  anode 
side.  The  electrode  height  is  used  as  characteristic  length  in  the 
calculation.  The  small  Peclet  number  values  imply  that  the  mass 
transport  rate  by  convection  is  much  lower  than  that  by  diffusion. 
Thus,  in  the  GDL  the  dominated  mechanism  of  mass  transfer  is 
diffusion. 


(a)  x-z  plane 


(b)  x-y  plane 

Fig.  5.  Velocity  fields  of  gas  mixture  in  the  anode:  (a)  x-z  plane  and  (b)  x-y 
plane. 

Fig.  7  shows  profiles  for  the  oxygen  mass  fraction  in  the 
cathode,  including  gas  channel  and  the  diffusion  layer.  At  low 
current  density  (Fig.  7a),  although  the  concentration  of  oxygen  is 
higher  in  the  channel  than  in  the  porous  medium,  oxygen  mass 
fraction  is  relatively  uniform  and  the  magnitudes  of  its  maxi¬ 
mum  and  minimum  are  in  the  same  order.  However,  the  oxygen 
mass  fraction  distribution  is  far  from  being  uniform  at  high  cur¬ 
rent  density  (Fig.  7b).  As  seen  from  the  figure,  the  minimum  of 
oxygen  concentration  is  located  at  the  corner  of  the  diffusion 
layer  of  outlet  and  its  magnitude  is  at  least  one  order  magni¬ 
tude  smaller  than  that  of  the  inlet  gas,  i.e.  the  maximum  mass 
fraction.  Such  oxygen  distribution  characteristics  imply  that  the 
local  current  density  is  non-uniformly  distributed  in  the  catalyst 
layer  at  high  average  current  density,  since  the  local  current  den¬ 
sity  is  dependent  on  the  oxygen  concentration  according  to  the 
Butler- Volmer  equation. 

Fig.  8  presents  the  distribution  of  local  current  density  in  the 
cathodic  catalyst  layer.  It  can  be  seen  that  the  distribution  is 
quite  uniform  at  low  current  density.  The  local  current  density 
is  somewhat  lower  over  the  shoulder  than  that  over  the  channel. 
Also  it  tends  to  rise  in  the  catalyst  layer  close  to  the  membrane 
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Fig.  7.  Oxygen  species  mass  fraction  (coQ,  non-dimension)  distribution  in  the 
cathode:  (a)  low  current  density  (7av  =  0.5  A  cm-2)  and  (b)  high  current  density 
(7av  =  1.2  A  cm-2). 


Fig.  6.  Velocity  fields  of  gas  mixture  in  the  cathode:  (a)  x-z  plane  and  (b)  x-y 
plane. 

where  the  local  activation  overpotential  is  increased.  However, 
over  the  entire  thickness  of  the  catalyst  layer,  the  local  current 
density  changes  comparatively  small.  While  at  high  current  den¬ 
sity,  the  distribution  pattern  becomes  very  different  from  that 
mentioned  above.  The  local  current  density  over  the  shoulder 
decreases  noticeably  compared  with  that  generated  in  the  area 
exposed  to  the  gas  channel.  The  minimum  current  density  is 
located  at  the  corner  of  the  catalyst  layer  over  the  shoulder  adja¬ 
cent  to  outlet,  where  the  oxygen  concentration  is  the  minimum. 
Its  magnitude  is  one  order  less  than  that  of  the  maximum,  which 
is  located  at  the  area  exposed  to  the  channel  inlet.  Also  it  is  vis¬ 
ible  that  the  local  current  density  of  the  aforementioned  corner 
is  increased  at  beginning  and  then  deceased  along  the  direction 
close  to  the  membrane.  The  phenomena  are  resulted  from  the 
mass  transfer  limitation  of  oxygen.  The  local  current  is  depen¬ 
dent  on  the  oxygen  concentration  and  activation  overpotential 
according  to  the  Butler- Volmer  equation.  The  oxygen  concen¬ 
tration  decreases  along  the  z-direction  close  to  the  membrane  due 
to  the  consumption  and  mass  transport  resistance,  while  the  acti¬ 
vation  overpotential  increases  along  the  same  direction  because 


of  the  facility  of  reaction.  At  the  beginning,  the  activation  over¬ 
potential  controls  the  electrochemical  reaction  and  the  current 
density  increases.  When  approaching  the  membrane,  the  oxygen 
concentration  becomes  so  small  that  the  reaction  is  controlled 
by  the  mass  transport  limitation.  The  local  current  density  then 
decreases  along  the  direction  close  to  the  membrane. 

Fig.  9  shows  the  contour  plots  of  cathodic  activation  overpo¬ 
tential  in  the  catalyst  layer  at  two  different  current  densities.  It  is 
seen  that  the  activation  overpotential  is  low  at  low  current  density 
and  more  or  less  uniformly  distributed  in  the  y-direction  over  the 
entire  cathode  catalyst  layer,  although  there  is  small  difference 
due  to  the  concentration  polarization  at  different  place.  At  higher 
current  density,  the  cathode  overpotential  is  higher  and  the  distri¬ 
bution  is  severely  non-uniform  in  the  y-direction,  influenced  by 
the  mass  transport  resistance  of  oxygen.  The  activation  overpo¬ 
tential  is  larger  over  the  shoulder  than  that  in  the  area  exposed  to 
the  channel.  Also  it  increases  downstream  along  the  channel  in 
the  direction  from  inlet  to  outlet.  This  phenomenon  is  apparent 
by  the  reason  of  the  poor  diffusion  of  oxygen  through  the  GDL 
and  the  decrease  of  oxygen  concentration  due  to  the  depletion  of 
reactants.  It  should  be  noted  that  the  local  activation  overpoten¬ 
tial  was  predicted  in  this  study.  However,  this  local  overpotential 
distribution  prediction  was  seldom  clearly  provided  in  the  exist- 
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(a)  Low  current  density  (lav=  0.5  A  cm 2) 
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(b)  High  current  density  (lav=  1 .2  A  cm  2) 


(b)  High  current  density  (lav=  1 .2  A  cm 2) 


Fig.  8.  Local  current  density  (7C,  Am-3)  distribution  in  the  cathode  catalyst 
layer:  (a)  low  current  density  (/av  =  0.5  A  cm-2)  and  (b)  high  current  density 
(/av  =  1.2  A  cm-2). 


Fig.  10.  Local  effectiveness  factor  (6,  non-dimension)  distribution  in  the  cathode 
catalyst  layer:  (a)  low  current  density  (/av  =  0.5  A  cm-2)  and  (b)  high  current 
density  (7av  =  1 .2  A  cm-2). 
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(a)  Low  current  density  (lav=  0.5  A  cm 2) 


(b)  High  current  density  (lav=  1 .2  A  cm 2) 


Fig.  9.  Local  activation  overpotential  (rjc,  V)  distribution  in  the  cathode  catalyst 
layer:  (a)  low  current  density  (7av  =  0.5  A  cm-2)  and  (b)  high  current  density 
(4V  =  1.2  A  cm-2). 


in g  modeling  papers,  such  as  [1-16].  This  feature  makes  our 
paper  more  clear  and  open  for  what  are  assumed  (or  adopted) 
and  what  are  predicted  during  the  calculation  of  the  overpoten- 
tial. 

The  local  effectiveness  factor  in  the  catalyst  layer  can  be 
predicted  in  this  simulation,  which  is  another  feature  of  this 
work.  Fig.  10  displays  the  effectiveness  factor  in  the  cathodic 
catalyst  layer  at  two  different  average  current  densities.  It  is 
seen  that  the  factor  is  less  for  the  current  density  of  1 .2  A  cm-2 
than  that  of  0.5  A  cm-2.  Also  the  factor  is  lower  in  the  region 
over  shoulder  than  that  in  the  region  exposed  to  channel,  which 
implies  that  the  diffusion  resistance  through  the  catalyst  particle 
is  greater  in  the  region  over  shoulder  than  that  in  the  region  over 
channel.  This  situation  is  further  intensified  at  higher  current 
density.  So  the  mass  transport  limitation  at  higher  current  density 
is  discovered  from  another  point  of  view. 

5.  Conclusion 

In  this  work,  an  isothermal,  steady-state,  three-dimensional 
multicomponent  computational  model  is  developed  for  PEM 
fuel  cell  with  straight  channels.  One  feature  of  the  model  is 
the  completeness  of  its  computational  domain  including  the  two 
flow  channels,  two  gas  diffusion  layers,  the  membrane  and  two 
catalyst  layers  with  their  physical  sizes.  Another  feature  of  the 
model  is  the  way  to  obtain  the  activation  overpotential.  The  local 
activation  overpotential  is  locally  predictedd  by  solving  the  elec¬ 
tric  potential  equations  of  membrane  phase  and  solid  phase,  not 
assumed  as  a  constant  through  the  catalyst  layer.  The  third  fea- 
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ture  is  that  the  agglomerate  model  is  introduced  into  the  3D 
model  to  account  for  the  effect  of  diffusion  resistance  through 
the  catalyst  particle,  and  the  local  effectiveness  factor  in  the  cat¬ 
alyst  layer  can  be  predicted. 

A  comprehensive  set  of  3D  continuity  equations,  momen¬ 
tum  equations  and  species  conservation  equations  and  potential 
equations  are  given  to  describe  the  flow  and  species  transport 
of  the  gaseous  mixture  and  electrochemical  kinetics  in  the  cou¬ 
pled  gas  channel  and  the  electrode.  They  are  solved  iteratively 
together  with  the  boundary  conditions. 

The  model  is  validated  by  comparing  the  polarization  curve 
of  model  prediction  with  experimental  data  of  Ticianelli  et  al. 
The  numerical  results  provide  the  detail  3D  velocity  vector  field 
and  pressure  field  of  gas  mixture.  Also  shown  are  the  3D  oxy¬ 
gen  concentration,  local  current  density  and  cathode  activation 
overpotential  distributions  at  different  current  densities.  The  dis¬ 
tribution  patterns  are  relatively  uniform  at  low  average  current 
densities  and  are  severely  non-uniform  at  higher  current  den¬ 
sity  due  to  the  mass  transfer  limitation,  which  implies  that  3D 
analysis  is  very  helpful  for  a  better  understanding  of  the  com¬ 
plicated  physical  process  in  a  fuel  cell.  From  the  predicted  local 
effectiveness  factor  in  the  catalyst  layer,  the  mass  transport  limi¬ 
tation  at  higher  current  density  is  indicated  from  another  point  of 
view. 

The  analysis  of  the  effects  of  operating  and  structural  param¬ 
eters  of  fuel  cell  on  its  performance  is  currently  in  progress  by 
using  this  model,  and  the  results  will  be  reported  elsewhere. 
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